function bspline_error

%%  Plot the max error at t=tmax 
%      as a function of the number of time points used
%  for the heat equation problem
%       diff(u,x,x) = diff(u,t)   for xL < x < xr, 0 < t < tmax
%  where
%      u = 0  at x=xL,xR  and  u = step at 0.25 & 0.75 at t = 0

clear *
clf

% get(gcf)
%set(gcf,'Position', [932 726 595 335]);
plotsize(595, 335)

% set parameters
tmax=0.1;
xL=0;
xR=1;

%%%%%%   plots as a function of time points  %%%%%%%%%%%%%

% iteration to determine the error as number of time points used is increased
N=20;

% generate the points along the x-axis, x(1)=xL and x(N+2)=xR
x=linspace(xL,xR,N+2);
h=x(2)-x(1);
hh=h*h;

ii=0;  M=2;
for i=1:11
	if i<8
		M=M+2;
	else
		M=2*M;
	end;
	ii=ii+1; points1(ii)=M;
	t=linspace(0,tmax,M);
	k=t(2)-t(1);
	lamda=k/h^2;
	errorC1(ii)=errorC(x,t,N,M,tmax,lamda);
end;	

ii=0;  M=2;
for i=1:14
	if i<12
		M=M+2;
	else
		M=2*M;
	end;
	ii=ii+1; points2(ii)=M;
	t=linspace(0,tmax,M);
	k=t(2)-t(1);
	lamda=k/h^2;
	errorBS1(ii)=bspline2(x,N,M,tmax,h,lamda);
end;	

% iteration to determine the error as number of time points used is increased
N=40;

% generate the points along the x-axis, x(1)=xL and x(N)=xR
x=linspace(xL,xR,N+2);
h=x(2)-x(1);
hh=h*h;

ii=0;  M=2;
for i=1:11
	if i<5
		M=M+2;
	elseif i<8
		M=M+3;
	elseif i<9
		M=M+6;
	elseif i<10
		M=M+3;
	else
		M=2*M;
	end;
	ii=ii+1; points3(ii)=M;
	t=linspace(0,tmax,M);
	k=t(2)-t(1);
	lamda=k/h^2;
	errorC2(ii)=errorC(x,t,N,M,tmax,lamda);
end;

ii=0;  M=2;
for i=1:14
	if i<5
		M=M+2;
	elseif i<10
		M=M+4;
	elseif i<13
		M=M+5;
	else
		M=2*M;
	end;

	ii=ii+1; points4(ii)=M;
	t=linspace(0,tmax,M);
	k=t(2)-t(1);
	lamda=k/h^2;
	errorBS2(ii)=bspline2(x,N,M,tmax,h,lamda);
end;	


% plot results
%subplot(2,1,1)
loglog(points1,errorC1,'--bo','MarkerSize',7)
hold on

loglog(points2,errorBS1,'-bd','MarkerSize',7)
loglog(points3,errorC2,'--rs','MarkerSize',7)
loglog(points4,errorBS2,'-r*','MarkerSize',7)
grid on
set(gca,'MinorGridLineStyle','none')
xlabel('Number of Time Points','FontSize',14,'FontWeight','bold')
ylabel('Error','FontSize',14,'FontWeight','bold')
set(gca,'FontSize',14)
legend(' N = 20 (CN)',' N = 20 (Bspline)',' N = 40 (CN)',' N = 40 (Bspline)',3);
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');
hold off
say=['Comparison of errors when u(x,0) has jumps'];
title(say,'FontSize',14,'FontWeight','bold')

% subfunction f(x,t)
function q=f(x,t)
q=0;

% subfunction g(x)
function q=g(x)
q=0.5*(sign(x-0.25)-sign(x-0.75));

% calculate error for bspline method
function q=bspline2(x,N,M,tmax,h,lamda)
aa=bcoeffs(N,M,x,lamda);
U=bsum(M,aa,x,h);
for ir=1:N+2
	exact(ir)=0;
	for ii=1:50
		an=2*(cos(0.25*pi*ii)-cos(0.75*pi*ii))/(pi*ii);
		exact(ir)=exact(ir)+an*exp(-pi*pi*ii*ii*tmax)*sin(pi*ii*x(ir));
	end;
end;
q=norm(U-exact,inf);

% error calculation for c-n
function q=errorC(x,t,N,M,tmax,lamda)
h=x(2)-x(1);  k=t(2)-t(1);
ue=crank(x,t,N,M,h,k,lamda);
for ir=1:N+2
	exact(ir)=0;
	for ii=1:50
		an=2*(cos(0.25*pi*ii)-cos(0.75*pi*ii))/(pi*ii);
		exact(ir)=exact(ir)+an*exp(-pi*pi*ii*ii*tmax)*sin(pi*ii*x(ir));
	end;
end;
q=norm(ue(:,M)-exact',inf);


function aa=bcoeffs(N,M,x,lamda)
% determine coeffs for initial condition '
aa=zeros(N,M);
a=4*ones(1,N); b=ones(1,N); c=b; z=zeros(1,N);
for i=1:N
	z(i)=6*g(x(i+1));
end;
aa(:,1) = tridiag( a, b, c, z )';
% determine coeffs for other t values '
a=(4+6*lamda)*ones(1,N); b=(1-3*lamda)*ones(1,N); c=b; z=zeros(1,N);
for j=2:M
	z=zeros(1,N);
	for i=2:N-1
		z(i) = (1+3*lamda)*aa(i-1,j-1) + (4-6*lamda)*aa(i,j-1) + (1+3*lamda)*aa(i+1,j-1);
	end;
	z(1)=(4-6*lamda)*aa(1,j-1) + (1+3*lamda)*aa(2,j-1);
	z(N)=(1+3*lamda)*aa(N-1,j-1) + (4-6*lamda)*aa(N,j-1);
	aa(:,j) = tridiag( a, b, c, z )'; 
end;


function ub=bsum(it,aa,xp,h)
% sum the series '
np=length(xp);
N=size(aa,1);
ub=zeros(1,np);
for ix=2:np-1
	sum=0;
	for kk=1:N
		xk=h*kk;
		xx=(xp(ix)-xk)/h;
		sum=sum+aa(kk,it)*bspline(xx);
	end;
	ub(ix)=sum;
	Lend=-aa(1,it)*bspline((xp(ix)+h)/h);
	Rend=-aa(N,it)*bspline((xp(ix)-h*(N+2))/h);
	ub(ix)=ub(ix)+Lend+Rend;
end;

function y=bspline(x)
% Calculate the value of cubic B-spline at point x '
% This is centered at x=0 and normalized so y=0 for x <-2 and x>2  '
x=abs(x) ;
if x>2,
  y=0 ;
else
  if x>1,
    y=(2-x)^3/6 ;
  else
    y=2/3-x^2*(1-x/2) ;
  end ;
end ;

% c-n method
function UC=crank(x,t,N,M,h,k,lamda)
UC=zeros(N+2,M);
for i=1:N+2
	UC(i,1)=g(x(i));
end;
a=-2*(1+lamda)*ones(1,N+2); b=lamda*ones(1,N+2); c=b; z=zeros(1,N+2);
a(1)=1; c(1)=0; a(N+2)=1; b(N+2)=0;
for j=2:M
	z(1)=0;  z(N+2)=0;
	for i=2:N+1
		z(i)=-lamda*UC(i+1,j-1)-2*(1-lamda)*UC(i,j-1)-lamda*UC(i-1,j-1)+k*f(x(i),t(j))+k*f(x(i+1),t(j-1));
	end;
	UC(:,j) = tridiag( a, b, c, z )';
end;
	
% tridiagonal solver '
function y = tridiag( a, b, c, f )
N = length(f);
v = zeros(1,N);   
y = v;
w = a(1);
y(1) = f(1)/w;
for i=2:N
    v(i-1) = c(i-1)/w;
    w = a(i) - b(i)*v(i-1);
    y(i) = ( f(i) - b(i)*y(i-1) )/w;
end;
for j=N-1:-1:1
   y(j) = y(j) - v(j)*y(j+1);
end;

% subfunction plotsize
function plotsize(width,height)
siz=get(0,'ScreenSize');
bottom=max(siz(4)-height-95,1);
set(gcf,'Position', [2 bottom width height]);